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Understanding the rheological properties of soft biological tissue is a key issue for mechanical 
systems used in the healthcare field. We propose a simple empirical model using Fractional Dynamics 
and Exponential Nonlinearity (FDEN) to identify the rheological properties of soft biological tissue. 
The model is derived from detailed material measurements using samples isolated from porcine 
liver. We conducted dynamic viscoelastic and creep tests on liver samples using a rheometer. The 
experimental results indicated that biological tissue has specific properties: i) power law increases in 
storage elastic modulus and loss elastic modulus with the same slope; ii) power law gain decrease and 
constant phase delay in the frequency domain over two decades; iii) log-log scale linearity between 
time and strain relationships under constant force; and iv) linear and log scale linearity between 
strain and stress relationships. Our simple FDEN model uses only three dependent parameters and 
represents the specific properties of soft biological tissue. 


I. INTRODUCTION 
A. Background 

Understanding the rheology —the study of materials 
with both solid and fluid characteristics in which the re¬ 
sponse to strain under applied stress is evaluated— of 
biological tissues is a key issue for current research in the 
human healthcare field. Rheology is relevant to many 
technological applications, ranging from biological sci¬ 
ence (e.g. medicine, sports, biology, biomechanics) to 
engineering (e.g., robotics, mechatronics, material me¬ 
chanics, control theory, computational mechanics, infor¬ 
mation technology). Recently, the healthcare field has 
realized the benefit of using intelligent machines (such 
as robots) that can physically interact with human. As 
a byproduct, the physical information measured by the 
machines can also be used for cyber system construction 
(such as machine learning). Understanding the physi¬ 
cal phenomena underlying the mechanical properties of 
human tissue has a great impact on bio-science and engi¬ 
neering. This knowledge will lead to further development 
of machines and systems in the healthcare held. 

Modeling of soft tissue rheological properties is a core 
technology for developing various healthcare machines 
and systems to assist human activity. For example, a 
mathematical model of target objects (human, organ, 
tissue, etc.) is required for mechanical design, motion 
planning, information processing, and machine/system 
control. These research and development areas require 
fundamental equations that are limited to the essential 
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properties of the macroscopic behavior of the target mat¬ 
ter (i.e., micro-scale modeling is not necessary). In short, 
the development of fundamental macroscopic models of 
the properties of biological matter are a key research issue 
pertinent to healthcare machines and systems designed 
for humans, organs, and tissues. 

In spite of their scientihc and technological impor¬ 
tance, mainly because they are difficult to model, very 
little knowledge has been established regarding the rheo¬ 
logical properties of soft biological tissues. The problem 
is due to a lack of established methods for sensing, pa¬ 
rameterizing, and information processing of rheological 
properties of soft biological tissue. The rheological prop¬ 
erties of soft biological tissue cannot be directly modeled 
in the same manner as synthetic matter. In general, an 
ordinary Linear Differential Equation (LDE) is used to 
model the target object. In other words, the terms of 
the equation for rheological properties have been gen¬ 
erally modeled using both ’Linear’ and ’Integer order’ 
differential equations explicitly or implicitly. However, 
the properties of soft biological tissues are different from 
synthetic matter and the LDE model tends to be inaccu¬ 
rate for data derived from experiments using soft biolog¬ 
ical tissue. Therefore, we aimed to develop a simplified 
fundamental macroscopic model of the rheological prop¬ 
erties of soft biological tissue that provides a good fit for 
experimental data. 


B. Goal and Motivation 

The goal of this study is to establish a universal fun¬ 
damental model to represent the macroscopic rheologi¬ 
cal properties of soft biological tissue, as well as a mea¬ 
surement method for these properties. Many researchers 
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have reported that the rheological properties of soft bio¬ 
logical tissues have distinct properties in comparison to 
industrial synthetic materials, such as metals. For exam¬ 
ple, researchers reported that biological tissues have vis¬ 
coelastic properties. Researchers have also reported that 
soft biological tissues exhibit a very nonlinear relation¬ 
ship between strain and stress. In this article, we describe 
the complex viscoelastic and nonlinear properties of soft 
biological tissue as ’rheological properties’. The moti¬ 
vation of this study is to propose a ’simple model’ that 
accurately represents the specific rheological properties 
of soft biological tissue. The model should be strongly 
correlated with experimental data derived from actual 
biological tissue. A ’simple model’ means that the model 
should utilize the minimum number of parameters, yield¬ 
ing a mathematical equation that is easy to understand 
and implement. Use of a simple model is essential for 
robust identification and discrimination of tissues using 
stress and strain information. 


C. Related research 

Numerous studies have dealt with both the non¬ 
linearity or/and viscoelasticity of biological tissue [I]- 
[22]. An Ordinary Linear Differential Equation (LDE) 
is generally used to model viscoelastic properties (e.g. 
Voigt/Maxwell/Kelvin model). However, LDE models 
do not fit data from biological tissues well. Eurthermore, 
a large number of parameters are used in LDE to in¬ 
crease model accuracy (e.g. generalized Voigt/Maxwell 
model), and these models only represent linear relation¬ 
ships between stress and strain. Hyperelastic models (e.g. 
Ogden, Moonine-Revlin model) are generally used to rep¬ 
resent stress-strain nonlinearity, although the number of 
parameters in hyperelastic models also tends to be large. 
Moreover, these models are time-independent and does 
not represent dynamic (viscoelastic) properties. Mod¬ 
els that neglect viscoelasticity and/or nonlinearity result 
in variability and incongruous analysis of the rheological 
properties of soft biological tissue. Specifically, parame¬ 
ter robustness decreases because the parameters readily 
change due to the duration of stress application and/or 
the magnitude of strain. 

The equations in the some related works have 

dealt with both viscoelasticity and nonlinearity. How¬ 
ever, these models tend to become overly complex and 
involve an excess number of material parameters to rep¬ 
resent these properties. We believe that a preferred 
model should have a small number of parameters that 
are strongly correlated with the experimental data. Ex¬ 
isting models with numerous parameters —such as those 
combining hyperelastic models with viscoelastic models- 
are unsuitable for identifying model parameters. The use 
of a large number of parameters leads to a risk overfit¬ 
ting of the parameter identification and ill-posedness of 
inverse problems. The large parameter number also in¬ 
creases computational costs. A model based on a simple 


equation with few parameters that is highly correlated 
with experimental data from soft biological tissues does 
not currently exist. 

Therefore, we have conducted studies aimed at devel¬ 
oping a model with these characteristics [23113(1] . The 
model is derived from comprehensive material data ob¬ 
tained from in vitro measurements of porcine liver iia- 
I26j . The model was also validated using in vitro breast 
tissue (fibroglandur tissue, fat, muscle) [571HH] and par¬ 
tially evaluated in muscle tissue |55J[3n] ■ The model com¬ 
bines a fractional differential equation with a polynomial 
expression for stress-strain nonlinearity, which consists of 
four parameters |23fl28j . However, the two parameters 
in the model —both parameters representing nonlinear 
properties— correlate and interfere with one another. In 
addition, the parameter identification from the experi¬ 
mental data of these two parameters is complex; specifi¬ 
cally, global searching and optimization is necessary. 


D. Objectives 

The objective of this article is to propose a sim¬ 
ple model that represents the rheological properties — 
meaning, viscoelastic and nonlinear properties— of soft 
biological tissues. Specifically, we propose a simple 
model, using only three dependent parameters, incorpo¬ 
rating fractional dynamics and exponential nonlinearity 
to identify rheological properties. The advantage of our 
model is that it is strongly correlated with various ex¬ 
perimental data and uses a small number of parameters, 
thereby rendering it suitable for parameter identification 
and inverse analysis. Eigure[2 shows an overview of this 
article. The model is derived from detailed material mea¬ 
surements using actual biological tissue. Specifically, we 
used samples isolated from various porcine livers. We 
selected liver samples because liver is a relatively simple 
tissue with low anisotropy when compared with other 
biological organs and tissues. We used a rheometer to 
measure the liver samples, as the rheometer can dynami¬ 
cally control and measure stress and strain applied to the 
sample. We conducted a dynamic viscoelastic test and 
creep test to derive and evaluate the model. Individual 
differences between liver samples -physical properties of 
biological tissues differ between individual samples- were 
represented by the values of model parameters. 


II. MATERIALS AND METHODS 

In this section, we explain how we measured and mod¬ 
eled the rheological properties —nonlinear and viscoelas¬ 
tic properties— of the samples. First, we introduce our 
rheological model scheme. We then explain the study 
materials and measurement procedures. 
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FIG. 1. The overview of this article. The model is derived from detailed material measurements using actual biological tissue. 
Specifically, we used samples isolated from various porcine livers. We used a rheometer to measure the liver samples in which 
the rheometer can dynamically control and measure stress and strain applied to the sample. Rheometers are commonly used 
to measure the fundamental physical properties of food samples. For example, we conducted a dynamic viscoelastic test and 
the creep test with several stresses to derive and evaluate our empirical model. Individual differences between liver samples 
-physical properties of biological tissues differ between individual samples- were represented by the values of model parameters. 


A. Proposed model 


The rheological model in this study relies on experi¬ 
mental data obtained from biological tissues. We first de¬ 
note the model equations ( [la| to enhance the read¬ 
ability of this article. The proposed rheological model 
utilizes Fractional Dynamics and Exponential Nonlinear¬ 
ity (FDEN). The equations are as follows: 




^iGx) = f 


{a; < Xb} ■ 
{a; > Xb} ■ 


(la) 

(lb) 


where x is strain (torsional strain), / is stress (torsional 
stress) and t is time, as variables; a is a non-integer 
derivative order representing the index of viscoelastic¬ 
ity, G is linear viscoelastic stiffness, is nonlinear vis¬ 
coelastic stiffness (Gi is dependent parameter), and Xb is 
the boundary strain in which the characteristics change 
to nonlinearity, as the parameters of the model, e is 
Napier’s constant. Each parameter should fulfill the fol¬ 
lowing relationship concerning the connectivity between 
linear (lal and nonlinear © equations -the exponential 
curve (lb) is a tangent to the straight line([I^-. 


1 

Xb — , 



( 2 ) 


The total number of parameters in the model is three 
(a, G and Gn as representative parameters) according to 


the relationships in ([^. The details of the experimental 
methods and derivation process of the model from the 
experimental data are described in the next sections. 


B. Materials and conditions 

Figure shows the details of the measuring compo¬ 
nents. We used porcine liver in the present study be¬ 
cause porcine abdominal organs have properties similar 
to those of humans, for example, the porcine abdominal 
organs are widely used in laparoscopic surgery training 
for novice surgeons. We chose to measure the proper¬ 
ties of liver samples because we thought that liver would 
be relatively easy to model -liver consists of homoge¬ 
neous and isotropic tissue-. We used cryogenically pre¬ 
served liver samples (4C on ice) that were taken within 24 
hours post-mortem and that did not include membranes 
or large blood vessels. Specimens were not frozen at any 
time during the procedure. We used a rheometer (AR550 
or AR-G2; TA Instruments, New Castle, DE) to measure 
the stress loaded on the sample and sample strain. The 
shear stress rheometer was selected because the shear test 
must be independent with the change of cross-sectional 
area in the stress calculation. In addition, the effect of 
gravity could be disregarded. From these measurements, 
the conventional shear strain x and conventional shear 
stress / were calculated, respectively. The details of the 
calculation are described in the Appendix]^ The liver 
sample was cut into slices (diameter 20 mm, height 5 
mm) and placed on a measurement table. The samples 
were soaked in saline solution at 35C during each test. 
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Sandpaper (P80 grain size) was attached to the top plate 
and the measurement table to prevent sliding. 


C. Procedures 

1. Initializing procedures 

After the saline solution reached the target tempera¬ 
ture, the gap was zeroed to the surface of the saucer. The 
saline solution was stable and there was no reflux flow. 
Each tissue sample was placed on a measurement table, 
and its thickness (=gap) was determined. The sample 
thickness was defined as the distance between the surface 
of the saucer and the surface of the parallel plate (part of 
the measuring device) at the time that the normal stress 
resulting from contact between the parallel plate and the 
sample reached 0.1 N. To engage the sample and parallel 
plate, preloading for 3 minutes and unloading for 3 min¬ 
utes was performed thrice under a load shear stress of 750 
Pa. The following series of experiments were conducted 
for each sample, after the above initializing procedures. 


2. Dynamic viscoelastic test 


3. Creep test and nonlinear measurement 

Torsional creep test was performed after dynamic vis¬ 
coelastic test. The creep test, in which step responses 
to strain are observed under constant stress, was repeat¬ 
edly performed, applying several stresses on the sample. 
Time series of strain data were measured during each ex¬ 
periment. The shear stress load ranged from 25 to 750 
Pa, and the time series of strain data were recorded for 
180 seconds at each stress level. Each test was performed 
at intervals of 180 seconds. The load shear stress during 
each interval was 0 Pa. The reference strain was set to 
0 at each creep test to account for residual stress and 
strain. We ignored the data obtained from 0 to 1 s be¬ 
cause of vibrations during the early transient stage. The 
details of this area are presented in our previous article 
[551 US]. Data were collected from 64 liver samples. 

III. RESULTS AND MODELING 
A. Mechanical complex impedance 

Here, mechanical complex impedance is defined as fol¬ 
lows: 


Sine-wave stress from 0.1 to 10 rad/s, providing 3% 
strain amplitude, was applied to the sample. The gain 
and phase delay of each frequency were measured. Then, 
gain (from stress to strain), phase delay, and mechanical 
complex impedance of the sample in each frequency were 
calculated. As shown in following experimental results 
(Fig§, 3% (= 0.03) strain amplitude is in the range 
where liver tissue exhibits linear responses. The effect 
of mass (inertia) and shear viscosity from the external 
normal saline solution could be disregarded at frequen¬ 
cies less than 10 rad/s. Data were collected from 6 liver 
samples. 



Radius^ 


FIG. 2. The details of the measuring components. We used 
porcine liver as sample of this study. We used a rheometer to 
measure the stress loaded on the sample and sample strain. 
The liver sample was cut into slices and placed on a measure¬ 
ment table. The samples were soaked in a saline solution at 
35 C during each test. Sandpaper was attached to the top 
plate and the measurement table to prevent sliding, d and R 
are the length and radius of the cylinder equation (A-1) and 
(A-2). R was 20 [mm] and d was 5 [mm] in the experimental 
setup of this paper. 


G* = G' +jG”. 


( 3 ) 


where j is unit imaginary number, G* is complex me¬ 
chanical impedance, G’ is storage elastic modulus and 
G” is loss elastic modulus. 

Typical experimental results of a dynamic viscoelastic 
test -in this section, mechanical complex impedance- of 
a sample are shown in Figj^ All liver samples exhibit 
the same trend as the typical sample; data trends are 
the same, however, model fit data and parameters are 
different. 

Both the storage elastic modulus G ’ and the loss elastic 
modulus G” increased with the frequency lu. We found 
that both storage elastic modulus G’ and loss elastic 
modulus G” exhibit a power law form over two decades. 
Furthermore, there is linearity in the log-log diagram in 
the change of G’ and G”, while the slopes of G’ and G” in 
the log-log diagram are nearly identical (approximately 
1/8=0.125). 

The mechanical complex impedance of our model has 
the same characteristics as the experimental results, i.e. 
power law form of G’ and G”, and same slopes of G’ 
and G”. The expansion of the equation to describe the 
explanation of the above characteristics is as follows. Our 
model is represented as equation Q -the same equation 
as (la) is described for readability- because the dynamic 


viscoelastic tests were conducted on a linear range in the 
stress and strain relationship. 


di^ 


{Gx) = f. 


( 4 ) 
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Because equation Q takes the form of a frequency 
transfer function, the complex shear modulus G* can be 
expressed in terms of the Laplace operator s as follows: 


represents loss elastic modulus. These parameters have 
the following relationship: 




( 5 ) 


Equation ([6| derived from the mechanical complex 
impedance of^^. 


G* = G Uujy 


( 6 ) 


G = 


^G'o^ + Gf. 

(8a) 

= Gcos(—a). 

(8b) 

= Gsin(—a). 

^2 ^ 

(8c) 


Euation (9a) and (9b I were derived from (7a) and (7b) 


by log-log transformation. 


Where uj is the 
(7a) and ([7b| with 
parts of (|^. 


frequency. Equation (§ expands to 
separation of the real and imaginary 


G' = G'fjuj^. (7a) 

G" = G"oUJ°‘ (7b) 


Where Gg is a constant parameter that represents stor¬ 
age elastic modulus and Gg is a constant parameter that 



FIG. 3. The typical experimental result of a dynamic vis¬ 
coelastic test -in this figure, mechanical complex impedance- 
of a sample. Blue plot is experimental result of storage elas¬ 
tic modulus G'. Red plot is experimental result of loss elas¬ 
tic modulus G. Both the storage elastic modulus G' and the 
loss elastic modulus G” increased as the frequency increased. 
Both storage elastic modulus G' and loss elastic modulus G” 
exhibit a power law form over two decades. Furthermore, 
there is linearity in the log-log diagram in the change of G' 
and G”. In addition, the slopes of G' and G” in the log- 
log diagram are nearly identical (approximately 0.125= 1/8). 
The G' of our model is blue line. G” of our model is red 
line. The G' and G” of our model, which parameters fit the 
typical experimental results, indicate that our model and the 
experimental results are highly correlated. 


log G' = a log UJ + log Gg . (9a) 

log G" = a logo;-|-log G'g. (9b) 


Thus, our model equation represents the trend in the 
experimental results, i.e. linearity on log-log diagram. 

The parameters (G and a) of equations (9a) and (9b) 
were identified by fitting the experimental results for each 
sample. We used the Extended Kalman Eilter (EKE) al¬ 
gorithm to identify the parameters (ref. Appendix C 
in detail) because equations (9a) and (9b) are nonlinear 
simultaneous equations -both equations include parame¬ 
ters (G and a )—. The G’ and G” in our model, which fit 
typical experimental results, are shown in Fig. showing 
that our model and the experimental results are strongly 
correlated. It should be noted that the order of deriva¬ 
tive a is not an integer; it is approximately 0.125 (= 1/8). 
Our model also fits all liver samples well. The coefficient 
of determination between our model and the experi¬ 
mental data from the frequency series of G’ and G” in 
all samples is about 90%. Tables |T| and [IT] list the model 
accuracy evaluation and fundamental statistics of model 
parameters. 


B. Bode diagram 

Typical experimental results of dynamic viscoelastic 
test -in this section, gain diagram and phase diagram - 
of a sample are shown in Fig|4} Figure]^ is a plot of the 
same data presented in Figj^the only difference being 
the expression of dynamic viscoelastic tests data. Gain 
-multiplicative inverse of G*- decreased as frequency in¬ 
creased. We found that gain assumes a power law form 
over two decades, indicating that there is linearity be¬ 
tween frequency and gain in the log-log diagram. We 
also found that phase delay remained constant over two 
decades. 

The bode diagram of our model has the same char¬ 
acteristics as the experimental result, namely power law 
form of gain and constant phase delay. The expansion 
of the equation to describe the explanation of the above 
characteristics is as follows. The Laplace operator of the 
bode diagram is as follows: 
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X(s) _ 1 _ 1 

F{s) Gs°‘ G{juj)°^ 


( 10 ) 


The model equations of gain, (11a) and (lib), are si m- 
ply calculated from (10) as follows. Equation (lib) is 
derived from log-log transformation of (11a): 


Gain = 


Gu^y 


Goj° 


\og{Gain) = — log G — a log w. 


(lla) 

(llb) 


In addition, the mode l equation of phase delay is de¬ 
rived from (7a) and (7b). 


G’ 1 TT 

Phase = arg( —) =arg(^^^^^^) = --a. (12) 

Thus, our model equation represents the trend ob¬ 
served in the experimental results. 

We calculated the gain and phase of our model via 
identihcation of the parameter of mechanical complex 
impedance for each sample because parameters were 
same. The gain and phase results from our model, which 
the parameters fit to the typical experimental results, are 
shown in Fig|^ which shows that our model and the ex¬ 
perimental results are strongly correlated. Our model 



Frcqucnc}' [rad/s] 


FIG. 4. The typical experimental result of a dynamic 
viscoelastic test -Bode diagram (gain diagram and phase 
diagram)— of a sample. Blue plots in gain diagram and phase 
diagram show the experimental result. Gain -multiplicative 
inverse of G*- decreased as the frequency increased. Gain 
assumes a power law form over two decades. Indicating that 
there is linearity between frequency and gain in the log-log 
diagram. We also found that phase delay remained constant 
over two decades. The gain and phase results of our model 
are red lines. Our model, which the parameters fit to the 
typical experimental results shows that our model and the 
experimental results are highly correlated. 


fit all liver samples well. Tables and list the model 
accuracy evaluation and fundamental statistics of model 
parameters. 


C. Creep test (Step response) 


Typical example of experimental results of a creep test 
-the creep response obtained by assuming the input step- 
stress- is shown in Figj^a). The strain of liver samples 
increased over a time interval of 180 s. Figurej^b) shows 
a log-log diagram of the same data described in Figj^a). 
We found the time series data of the creep response exhib¬ 
ited a power law form over two decades. This indicates 
that there is linearity between time and strain in the log- 
log diagram (Figj^b)). A model equation of strain in the 
creep test can be calculated. We assumed that equation 
([^ is valid for a single creep test, while nonlinearity was 
evaluated by a series of creep tests under sev eral a pplied 
stresses. Specifically, equation 0 becomes ( |13a| ) if 0 
is solved for the conditions of the creep test. Here, the 


applied stress is constant fc- Equat ion (13b) is derived 
from log-log transformation of (13a). 


loga; = alogt-I-logXc- (13b) 

where x is strain and t is time; fc is constant stress; G 
is proportional factor representing stiffness at each strain, 
as parameter r() is the gamma function. Xc is the coeffi¬ 
cient determining the strain value as a parameter, which 
is defined as follows: 


Xc = 


fc 


GT(1 -b a) ■ 


(14) 


In this case, the Riemann-Liouville definition (15) -but 


not only this definition- was used to solve the fractional 
integration of 0. 


Z7-“/(t) = 


r(a) 


J 0 


/(OdC- (15) 


Thus, our model equation represents the trend ob¬ 
served in the experimental results. The parameters (a, 
Xc) of equation (13b) were identihed by fitting the ex¬ 
perimental results in the log-log domain. We used the 
LSM algorithm to identify the parameters of equation 
(13b) -linear regression- for each sample. We calculated 


the other independent parameter G via equation (14). 


The time series displacement data from our model, the 
parameters of which ht the typical experimental results, 
are shown in Figsj^a) and (b). Figsj^a) and (b) show 
that our model and the experimental results are strongly 
correlated. Our model fit all liver samples well. The co¬ 
efficient of determination between our model and the 
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experimental data from the time series of displacement 
in all samples at all stresses exceeded 99%.Tables |l] and 
[IT] list the model accuracy evaluation and fundamental 
statistics of model parameters. 


D. Nonlinearity measurement 

Nonlinear properties of samples were investigated 
based on a series of creep tests under several applied 
stresses. Specifically, we looked at the relationship be¬ 
tween the constant applied stress fc and the strain coef¬ 
ficient Xc in a series of creep tests using several stresses. 
Typical experimental results for nonlinearity measure¬ 
ment of the sample are shown in Figj^a) and (b). Figure 
[^b) shows a semi-log diagram of the same data described 


(a) 0 

O.I 

s 

1 0 
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FIG. 5. (a) A typical example of experimental results of 

a creep test -the creep response obtained by assuming the 
input step stress-, (b) shows a log-log diagram of the same 
data described in (a). The Blue plot in (a) and (b) are ex¬ 
perimental results. The strain of the sample increased over a 
time interval of 180 s. We found the time series data of the 
creep response exhibited a power law form over two decades. 
This indicates that there is linearity between time and strain 
in the log-log diagram. The time series data of our model are 
red lines. Our model equation represent the trend observed in 
the experimental results. These hgure shows our model and 
the experimental results are highly correlated. 


in Fig^a). Figure |^a) shows the relationship between 
Xc and fc exhibit linear characteristics under low strain 
conditions. The stress nonlinearly increased under high 
strain conditions. Figure [^b) shows that the stress in¬ 
crease during high displacement is linear in the semi log 
scale; stress increases exponentially in the linear scale 
space. We found that linear straight line and nonlinear 
curves connected smoothly, with the exponential curve 
tangent to the straight line region. We modeled nonlin¬ 
ear properties of the soft biological tissue based on t hese 
result s and consi derat ions, as shown in equation! 16a I and 
(16b I. Equation (16c I is derived from log-log transforma¬ 
tion of (|16b|). 


o 

II 

{Xc < Xb}. 

(16a) 

= fc 

{Xc > Xb}. 

(16b) 

GnXc + log Gi = log fc 

{Xc > Xb}. 

(16c) 


where Xc is strain and fc is stress. G is a proportion¬ 
ality factor (referred to as linear viscoelastic stiffness or 
linear stiffness herein), Xb is the boundary strain between 
the linear and nonlinear range (called boundary strain), 
Gn is a proportionality factor in the log space (called 
nonlinear viscoelastic stiffness or nonlinear stiffness), Gi 
is the dependent parameter. Each parameter should ful¬ 
fill the following relationship due to the exponential curve 
(16b) being a tangent to the straight line (|16a|. Details 


of this relationship are shown in Appendix |B 


Xb = 


G, = 


1 

G 

Gri.e 


(17a) 

(17b) 


The parameters {a, G, G„, Gi and Xb) of equations 


(16a I and (16c) were identified by fitting the experimen¬ 


tal re sults for e ach sa mple. Thu s, the equation (16a) and 
(|16c[) becomes (18a) and (18b) when the only indepen¬ 


dent parameters are used. 


Gxc = fc {xc<-p^}- (18a) 

G„a;c-blog(-^) =log/c {xc> ^}- (18b) 

We used the Extended Kalman Filter Algorithm to 
identify the parameters of equation (18a) and (18b) (ref. 
Appendix in detail), as they are nonlinear simulta¬ 
neous equations. The stress-strain relationship of our 
model, which fit the typical experimental results, are 
shown in Figj^ show that our model and the experimen¬ 
tal results are strongly correlated. Our model also fit 
all liver samples well. The coefficient of determination 
Ft? between our model and the experimental data in all 
samples was about 95%. Tables [I] and [T^ list the model 




















accuracy evaluation and fundamental statistics of model 
parameters. 


Thus, we derived the nonlinear equations for our model 
shown in (la) and (lb). We assume here that equations 
(16a I and (16b I can hold true more generally in the ab¬ 
sence of creep tests. 


(a) 



Siniiii 


FIG. 6. (a) shows the relationship between strain Xc and 

stress fc exhibit linear characteristics under low strain condi¬ 
tions. The stress nonlinearly increased during high strain con¬ 
ditions. The stress nonlinear ly increased during high strain 
conditions, (b) shows that the stress increase during high 
displacement is linear in the semi log scale; stress increases 
exponentially in the linear scale space. We found that linear 
straight line and nonlinear curves connected smoothly, with 
the exponential curve tangent to the linear straight region. 
The stress-strain relationship of our model, which fit the typ¬ 
ical experimental results, are shown in red line. Our model 
and the experimental results are highly correlated. There is 
an overlapping area where both linearity and log scale lin¬ 
earity. Linearity holds to a certain degree over the boundary 
strain Xb- log scale linearity holds to a certain degree before 
the boundary strain Xb- 


IV. DISCUSSION 


The main contribution of this article is the proposal 
of a Fractional Dynamics and Exponential Nonlinearity 
(FDEN) model to identify the rheological properties of 
soft biological tissue. We found from experimental results 
that biological tissues have specific properties: i) power 
law increases in storage elastic modulus and loss elastic 
modulus of the same slope; ii) power law gain and con¬ 
stant phase delay in frequency domain over two decades; 
iii) log-log scale linearity between time and strain rela¬ 
tionships over two decades; and iv) linearity in low strain 
range and log scale linearity in high strain range between 
strain and stress relationships. The FDEN model uses 
only three dependent parameters (such as a, G and G„) 
and represents the specific properties of soft biological 
tissues. The advantage of our model is that it strongly 
correlates with various experimental data, as shown in 


the section III In addition, the small number of parame¬ 
ters used is valuable because it is suitable for parameter 
identihcation and inverse analysis. For example, the pa¬ 
rameter identification methods in this article are basic, 
with only the Least Square Method (LSM) and Extended 
Kalman Eilter (EKF) being used. Lastly, the meaning of 
each parameter is intuitively understood (for example, a: 
ratio of viscoelasticity, G: linear stiffness, Gn- nonlinear 
stiffness) and it is possible to compare the values with 
other tissue. The details of the discussion are described 
in the following sections. 


A. Viscoelastic model using fractional calculus 

We found from experimental results that soft biolog¬ 
ical tissues have specific properties, as described above 
in i)-iv). Single terms in the fractional dynamics model 
0 represent the specific viscoelastic properties of soft 
biological tissues. Fractional calculus is an approach to 
mathematically describe natural phenomena that are re¬ 
lated to viscoelastic behavior [ST]. 

Fractional calculus is a branch of mathematical analy¬ 
sis concerned with taking real or complex number pow¬ 
ers of differential operators. Fractional dynamics is a 
field of study in physics and mechanics concerned with 
investigating the behavior of objects and systems that 
are characterized by power law non-locality, power law 
long-term memory or fractal type properties by using in¬ 
tegration and differentiation of non-integer orders, i.e., 
by fractional calculus methods [35] . Fractional dynamics 
models have proven to be powerful tools in describing the 
dynamic behavior of various materials. The advantages 
of fractional dynamics models are their ability to de¬ 
scribe real dynamic behavior and the fact they are simple 
enough for engineering calculations [33] . The equations 
for rheological models are generally based on stress-strain 
analyses and are traditionally represented with deriva¬ 
tives of integer order (ordinary differential equation). In 
other words, traditional methods to fit the viscoelastic re- 
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sponse include several spring and dashpot elements. Re¬ 
cently, fractional dynamics models proved to be efficient 
in describing rheological materials such as rubber and 
tissues, reducing the number of parameters and showing 
a power law response [M] . 

Over the last years, fractional calculus has become an 
important tool in the analysis of viscoelastic materials 
composed of synthetic polymers [3S]. For example, Ca- 
puto et al. found good agreement with experi¬ 

mental results when using fractional calculation for the 
description of viscoelastic materials and established the 
connection between fractional calculation and the theory 
of linear viscoelasticity [35] ■ Several authors [301111] have 
also suggested the use of differential or integral equations 
of fractional order to describe viscoelastic behavior that 
is intermediate between purely elastic and purely viscous 

m- 

Although fractional calculus has wide application in 
describing the solid-liquid duality of synthetic polymers, 
it had until recently attracted limited attention in the 
field of biological materials, biomechanics, bio-rheology 
and cell viscoelasticity. Suki et al. [33] found the pres¬ 
sure/volume response of a whole lung to be characterized 
by fractional calculus. Fractional calculus is also use¬ 
ful in bio-fields because many tissue-like materials (poly¬ 
mers, gels, emulsions, composites, and suspensions) ex¬ 
hibit power law responses to applied stress or strain [33] ■ 
Although such models are widely used for synthetic ma¬ 
terials, they have been progressively substituted in the 
community by models relying on fractional calculus for 
biological materials [33]. Yuan et al. [31113, studied 
lung tissue and found its fractional order of evolution, 
while Chen et al. [31] applied the same model to agarose 
gels used for culturing tissues, particularly cartilage. An 
example of the power law behavior of elastic tissue was 
observed recently for viscoelastic measurements of blood 
vessels, where the analysis of these data was most con¬ 
veniently performed using fractional order viscoelastic 
models [33] ■ Recently, the framework of fractional calcu¬ 
lus has also been used in the research of magnetic reso¬ 
nance elastography [33] |3S] ■ As above, fractional dynam¬ 
ics are gaining popularity in the field of viscoelasticity, 
with data and models already reported for the liver [H- 
[26], breast tissues [28], lung [331139], vessels [33], muscle 
[351 [35] . muscle cells [30], tendons [33 and blood cells 
[53] . In short, research on fractional calculus has been 
applied widely to many fields, including biological mate¬ 
rials. 

The parameter a in fractional equation @ is in the 
order of a derivative that is commonly taken to range 
between 0 and 1. If a is 0, equation Q describes the 
behavior of a spring where G specifies the springs stiff¬ 
ness. If a is I, equation defines a dashpot, in which 
G defines the viscosity. Thus, the fractional equation Q 
interpolates between the material behavior of a spring 
and that of a dashpot [35]. The rheological element that 
refers to equation 0 was therefore introduced by Koeller 
and termed a ’springpot’ [331133] ■ As such, the derivative 


order a represents the index of viscosity of the system in 
the fractional dynamics model. A viscoelastic material is 
more governed by elastic properties than by the viscous 
properties when the derivative order a is close to 0. A 
viscoelastic material is more governed by viscous proper¬ 
ties than by elastic properties when the derivative order 
a is close to I. 

The value of the derivative order a was 0.125 (= 1/8) 
from the experimental results displayed in this article, in¬ 
dicating that the characteristics of soft biological tissue 
(liver) are intermediate between those of elastic and vis¬ 
cous bodies and is relatively close to an elastic material. 


B. Fractional calculus for the dynamic viscoelastic 
and creep tests 


In this article, the viscoelastic properties of soft biolog¬ 
ical tissues (liver) have been examined. The simple em- 
piric al equ ations describing strain creep (equation( ]13a| ) 
and (13bl) have been put in a concise mathematical 
framework. We have chosen to describe viscoelasticity in 
terms of fractional calculation as in equation Q . Certain 
important advantages must be emphasized, namely, frac¬ 
tional calculus has the following advantages [331133] ; i) 
fractional dynamics models accurately describe complex 
models with fewer parameters, ii) they improve curve fit¬ 
ting, principally with power law responses, iii) they al¬ 
lowed a physical justification in rheological and tissue 
samples. 


1. Dynamic viscoelastic test 

Measurements of mechanical complex impedance G* 
in dynamic viscoelastic tests over a wide range of forcing 
frequencies (I0“^-I0^ rad/s) in tissue sample revealed 
that the frequency dependence of rheological behavior 
represents a weak power law relationship over a wide 
range of frequencies —For example, Fabry |54j reported 
that a weak power law relationship held over a range of 
frequencies (I0“^-I0^ Hz) in muscle cells—. The storage 
modulus G’ increases with increasing frequency accord¬ 
ing to a weak power law with a power law exponent of ap¬ 
proximately 0.125. The loss modulus G” also follows the 
power law with a power law exponent of approximately 
0.125. Fractional calculus provides a natural framework 
for describing such weak power law relationships E3- 

In contrast, mechanical models using ordinary differ¬ 
ential equations have long been used, and their qualita¬ 
tive behavior is not representative of the actual behavior 
of materials. The characteristics of the frequency de¬ 
pendencies could be similar; however, the slopes of the 
experimental results do not fit those of the theoretical 
curves [33]. The shortcomings of ordinary differential 
models can be recognized by comparing the frequency 
curves exhibited by a material with those predicted by 
the models. The weak power law behavior cannot be ac- 
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counted for by standard viscoelastic models characterized 
by ordinary differential equations [51] : i) storage modu¬ 
lus G’ should remain constant at low frequencies, which 
would indicate elastic behavior in ordinary differential 
equations and ii) loss storage modulus G” increases and 
approaches a power law exponent of 1 at high frequen¬ 
cies, indicating viscous behavior in an ordinary differen¬ 
tial equation. 


2. Creep test 

Recent studies also indicated that the time domain 
data of tissues are well represented by a simple empirical 
equation involving a power law in time. Some studies also 
reported that creep responses represent power law stress 
to a step input in the time domain. Fung [1] demon¬ 
strated in his theory that a distribution of time constants 
proportional to power of time over a finite range of time 
constants is appropriate for many tissues |4!1| . Djordjevic 
and co-workers [51] reported that a parallel combination 
of a fractional calculus (springpot) and a dashpot prop¬ 
erly predict the measured values for a rheological model 
of cultured smooth muscle cells [U]. As above, fractional 
calculus provides a natural framework for describing such 
power laws in the time domain [51] . 

In contrast, mechanical models using ordinary differ¬ 
ential equations lack consistency between their qualita¬ 
tive behavior and the real behavior of materials curves. 
Although the characteristics of time dependences could 
be similar, the slopes of experimental results do not fit 
those of the theoretical curves. The shortcomings of or¬ 
dinary differential models can be recognized by compar¬ 
ing the time domain curves observed for a material with 
those predicted by the models [33]. The power law be¬ 
havior cannot be accounted for by standard viscoelas¬ 
tic models characterized using ordinary differential equa¬ 
tions, such that; i) strain should remain constant at suffi¬ 
cient elapsed times, which would indicate elastic behavior 
in ordinary differential equations and ii) exponential in¬ 
creases in transient state, which would indicate viscous 
behavior in ordinary differential equations. 


C. Fractal structure and the Fractional ladder 
model 

Theoretical aspects of the fractional structure of soft 
biological tissue are partially explained with tissue frac¬ 
tal geometry in nature and its relationship to fractional 
calculus. Currently, fractal geometry and fractional cal¬ 
culus are applied to phenomenological theories for com¬ 
plex systems m- Soft biological tissues also have fractal 
structures, such as in Figj^ A fractal is a natural phe¬ 
nomenon or a mathematical set that exhibits a repeating 
pattern displayed at every scale. For example, Schiessel 
and Blumen [35], Schiessel et al. [55] and Heymans and 
Bauwens [23] have demonstrated that fractional equa¬ 


tions, such as Q, can be realized physically through the 
fractal structure of hierarchical arrangements of springs 
and dashpots like ladders [55]. A related work [57] de¬ 
scribes the details of the following layered fractal models 
of soft biological tissues based on the schema shown in 
Figj^ The first panel (a) displays an infinite number of 
thin elastic membranes and viscous compartments. The 
second panel (b) is a magnified region of the first panel 
(a) and the third panel (c) is a magnified region of the 
second panel(b), showing the self-similar layered struc¬ 
ture. By allowing the number of structural components 
to extend indefinitely, the self-similarity of biological me¬ 
dia is revealed. This topology is also depicted in Fig. [^ 
where the alternating elastic and viscous components are 
visualized as a self-similar hexagonal packing of spheres 
within spheres. In order to capture these fractal compo¬ 
nents with the elastic membranes and viscous saline of bi¬ 
ological tissue, a fractal ladder of springs and dashpots is 
described in Fig Em- A paper m described the prop¬ 
erties of Fig. which presents the fractional deriva¬ 

tive term with the order of derivative 1/2. Similar fractal 
tree networks were considered m to model the other or¬ 
der of fractional calculus (orders of 1/4 and 1/8 are also 
described in Fig. |^b) and (c), respectively). These re¬ 
cursive ladder expansions provide various parameters of 
derivative order. Specifically, the ladder model can also 
be considered as a fundamental mechanical component of 
fractional derivative term, allowing more complex fractal 
networks, or recursive ladders, to be constructed. For in¬ 
stance, consider a recursive ladder model constructed by 
replacing the viscous damper in Fig. |^a) with a fractal 
ladder, producing the arrangement shown in Fig. ib) 
with the order of derivative 1/4. Similarly, a recursive 
ladder may be constructed by replacing the springs in 
Figj^b) with a fractal ladder, producing the arrange¬ 
ment shown in Figj^c) with the order of derivative 1/8. 

The parameter a (order of derivative and also power 
law exponent) was about 1/8 (= 0.125) according to the 
experimental results of the dynamic viscoelastic tests in 
Figs.([^ and Q. This suggests that liver tissue has a 
complex fractal structure, such as in Figj^c), where the 
fractional ladders are several times renormalized. 


D. Main contribution 

From the viewpoint of fractional calculus, the main 
contribution of this study is to propose integration with 
fractional calculus and nonlinear equations, in other 
words, utilizing the nonlinearity of fractional calculus to 
describe soft biological tissue. The history of fractional 
dynamics in viscoelastic materials is long standing; how¬ 
ever, the nonlinearity of the fractional equation for soft 
biological tissue has not yet been considered in related 
studies. The theoretical contribution of the present study 
is to empirically reveal the order of fractal structure of 
soft biological tissues (Figj^c)) via the analysis of the 
derivative order {a= 1/8). 
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We believe that the fractional model with a single term 
is most suitable for the identification of bio-rheological 
properties, while many previous studies also have pro¬ 
posed serial and/or parallel arrangements of ordinary or¬ 
der models and fractional order models (such as a frac¬ 
tional generalized Voight model with a fractional term). 
The single fractional derivative term has the strong ad¬ 
vantage of high model accuracy - although the frequency 
range was relatively low in this study- and the power 
law relationship is suitable for parameter identification, 
as described in the following section[TVF[|IVG]and[lVH| 


E. Exponential nonlinearity 


Figure|^and equation (16b I show that stress on soft bi¬ 
ological tissue increases exponentially with strain. These 
trends are generally known in fields such as economics 
and in natural evolutionary processes. For example, 
value grows exponentially with time, technology has ad¬ 
vanced at an exponential rate |58j (exponential growth 
of computing power is known as Moores law), market 
price in inflation shows exponential growth [59] and pop¬ 
ulation growth (such as Malthusian Theory of Popula¬ 
tion) is exponential. The experimental results of stress 
and exponential models imply that the behaviors in the 
stress-strain relationship may have a similar mathemat¬ 
ical structure, although the variable is not time, but, 
strain. In this theory, the exponential growth evolves due 
to a linear positive feedback mechanism, such as equation 
(19b); an upward change in stress induces further stress 


Extracellular Matrix Cytoplasm Nucleoplasm 



FIG. 7. Fractal structure of soft biological tissues [57|. A 
repeating pattern with thin elastic membranes and viscous 
components is displayed at every scale. The first panel dis¬ 
plays an infinite number of thin elastic membranes and viscous 
compartments. The second panel zooms in on the first panel, 
thus showing the self-similar layered structure. By allowing 
the number of structural components to extend indehnitely, 
the self-similarity of biological media is revealed. 


increases rather than just incremental additons. 



{x > Xb] A {cc < —Xb] 


(19a) 

(19b) 


Where is the intermediate variable between the up¬ 
per and lower equations, which is related to force. We in¬ 
troduce a second-order partial differential equation here 
because of the negative and positive symmetry properties 
of the strain and stress relationship, as shown in Figure 
1^ and equation (201. The solution for the second-order 
order partial differential equation is as follows; 


GiC = / {a: < -Xb}. 

(20a) 

Gx = f {—Xb < X < Xb}. 

(20b) 

= f {x>xb}. 

(20c) 


The solution for the first-order partial differential equa¬ 
tion only represents positive or negative exponential 
stress changes. 

The smoothness of a fundamental mathematical func¬ 
tion largely affects the robustness of identification, in¬ 
verse analysis, structure analysis in computer simula¬ 
tions, and optimization using structure analysis. In par¬ 
ticular, the smoothness at a point where mathematical 
function changes (specifically, x = Xb) is important. Dif¬ 
ferentiability class is generally used in the classification 
of functions according to smoothness, more specifically, 
the properties of their derivatives. The differentiability 
class of our nonlinear model - between equations such 
as (16a) and (16b), (20a) and (20b), (20b) and (20c)- 
is Cl a.t X = Xb,—Xb] this means curves are continuous 
and differentiable at x = Xb- In other words, curves are 
joined and first derivatives are continuous at a; = Xb, —Xb- 
Thus, the smoothness of function in our nonlinear model 
is maintained when compared with linear models. 

These smooth characteristics are an advantage of our 
model because the robustness of calculations in the 
boundary between linear and nonlinear characteristics is 
high. We can estimate the origin of strain from data in 
the nonlinear range due to the constraint condition of 
parameters described in equation (17a). Equation (17a) 
refers to subtangent (xb) as a word in geometry is con¬ 
stant in exponential function. The subtangent can be 
estimated using the tangent {Gn) of the semi-log graph 
of the stress-strain relationship. This is an important 
characteristic because the zero strain point of soft bio¬ 
logical tissue is generally difficult to define. The zero 
strain point is difficult to define because: 1) the zero 
strain point cannot be defined from linear data; 2) defor¬ 
mation of soft biological tissue is relatively large and is 
markedly affected by gravity force; and 3) viscoelasticity 
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FIG. 8. Fractal ladder of springs and dashpots in order to capture these fractal components with the elastic membranes and 
viscous saline of biological tissue. Fractal tree networks were considered to model fractional calculus (orders of 1/2, 1/4 and 
1/8 are described in (a), (b) and (c), respectively). These recursive ladder expansions provide various parameters of derivative 
order. Specifically, the ladder model can also be considered as a fundamental mechanical component of fractional derivative 
term, allowing more complicated fractal networks, or recursive ladders, to be constructed. The parameter a (order of derivative 
and also power law exponent) was about 1/8 (= 0.125) according to the experimental results of the dynamic viscoelastic tests 
in Figs.Q and This result suggests that liver tissue has a complex fractal structure such as in (c), where the fractional 
ladders are several times renormalized. 


of soft biological tissue results in difficulties measuring 
the zero strain point. 



FIG. 9. Negative and positive symmetry property of strain 
and stress relationship to introduce second order partial dif¬ 
ferential equation(in the case of G = 1000 ad Gn = 10). The 
solution of first order partial differential equation only rep¬ 
resent positive or negative exponential stress changes. The 
solution of second order partial differential equation repre¬ 
sent positive and negative exponential stress changes such as 
this figure. 


F. Time and frequency scale invariant 

One attribute of power laws is their scale invariance. 
Scaling the argument by a constant factor causes only 
a proportionate scaling of the function itself. Scaling 
by a constant simply multiplies the original power law 
relationship by the constant (parameter a in this article). 
Thus, it follows that all power laws with a particular 
scaling exponent are equivalent up to constant factors, 
as each is simply a scaled version of the others. There 
is no internal time scale that could typify the dynamics 
and no time characteristics are evident. 

Time scale invariance during creep tests can be de¬ 
scribed based on the above discussion, and creep response 
is not tied to any time scale; thus, it may be regarded 
as being scale-free. More specifically, as explained in this 
article, obtaining numerous time series data from creep 
tests is not necessary due to invariance in the time scale 
property. Only two data points at any time point are 
sufficient to robustly identify the parameter of equation 
Naturally, this is only a theory pertaining to iden¬ 
tical conditions, and many data points are preferable to 
enhance the robustness of measurements. 

Frequency scale invariance during dynamic viscoelastic 
tests can be described in the same manner. There are no 
internal frequency scales that could typify the dynamics 
and no characteristic frequency was evident. Mechani¬ 
cal impedance responses are not tied to any frequency 
scale; thus, they may be regarded as being scale-free. 
































































13 


More specifically, obtaining numerous frequency series 
data from dynamic viscoelastic tests are not necessary 
due to invariance in frequency scales. Only two data 
points at any frequency point are sufficient to robustly 
identify the parameter of equation Q. Naturally, this 
is only a theory pertaining to identical conditions, and 
many data points are preferable to enhance the robust¬ 
ness of measurements. 


G. Strain scale invariance 

Strain scale invariance -while nonlinearity is not power 
law- also holds at the linear area and nonlinear scale. 
The relationship between stress and strain in the linear 
range exhibits strain scale invariance because of linear¬ 
ity. The relationship between the logarithms of stress and 
strain in the nonlinear range also exhibits strain scale in¬ 
variance because of the linearity of the semi-log space. 
This scale invariance produces a strong relationship for 
identifying parameters. Theoretically speaking, a two 
point data set in the linear range (including stress = 0 
and strain =0) is sufficient to identify linear stiffness G 
(slope of stress and strain in the linear space). In addi¬ 
tion, a two-point data set in the nonlinear range is suf¬ 
ficient to identify nonlinear stiffness G„ (slope of stress 
and strain in the semi-log space). Of course, the afore¬ 
mentioned data set numbers are theoretical for identical 
situations, many data points are preferable to enhance 
measurement robustness. 

In addition, strict classifications between linear and 
nonlinear range should not be neccesary because of the 
smooth connectivity between the linear and nonlinear 
properties of soft biological tissues (Fig. [^. There is an 
overlapping area where both linearity and log scale linear¬ 
ity are held (about strain from 0.10 to 0.15 in Figj^. Lin¬ 
earity holds to a certain degree over the boundary strain 
Xb- Furthermore, log scale linearity holds to a certain 
degree before the boundary strain Xb- In other words, 
the data set in the overlapped area includes information 
from both the linear and nonlinear range. These proper¬ 
ties, which are common to soft biological tissues and the 
FDEN model, are useful for identifying parameters and 
it is possible to identify both linear and nonlinear stiff¬ 
ness using only data sets derived from the overlapping 
area. However, further research is required to confirm 
this hypothesis. 


H. Identification algorithm 

Acceding to the simple model equation, its power 
law properties and its semi-log scale linearity, parame¬ 
ter identification in the FDEN model is simpler when 
compared to other reported models. The small number 
of parameters in the FDEN model contributes to a simple 
algorithm and parameter identiheation. These character¬ 
istics may be also effective in inverse analysis of computer 


structural simulations of tissue/organ deformation. For 
example, the parameter identification methods in this ar¬ 
ticle are basic, with only the Least Square Method (LSM) 
and Extended Kalman Filter (EKF) being used. The pa¬ 
rameters of the model in the creep test can be identified 
using LSM. EKF was used to identify the parameter of 
dynamic viscoelastic tests -mechanical impedance- be¬ 
cause equations (9a I and (9b I are nonlinear simultaneous 
equations. EKF was also used to identify the parameter 
of nonlinear models between stress and strain relation¬ 
ship because equations (I8a) and (18b I are nonlinear si¬ 
multaneous equations. 

The parameters of the model in the Bode diagram 
were identified via parameter identification of mechanical 
impedance in this article. It should be noted that LSM is 
sufficient to identify the value of a and G in (11b). The 
correlation between model (121 and experimental data 
of the phase should be carefully checked in the case, as 
phase data can affect identiheation of the parameters. 


I. Parameter variation 

Tables |T] and [IT] list the model accuracy evaluation and 
fundamental statistics of model parameters. Statistical 
analysis of the samples used in this study revealed that 
the maximum values of a and G and were approxi¬ 
mately 1, 4 and 3 times the minimum values, respectively. 
These results indicate that the linear stiffness G and the 
nonlinear stiffness G„ (also, Xb as an independent pa¬ 
rameter) has a large degree of variation when compared 
with ratio of viscoelasticity a. The one of limitation of 
this study is that the number of samples was insufficient 
to stastically analyze the variations in each parameter. 
We plan to study other tissue types in order to compare 
and discriminate between tissues using these static and 
variable parameters (for example, |28]L 


J. Limitations 

The main limitation of this study is that we only mea¬ 
sured and evaluated liver samples. Similar evaluations 
must be performed with other tissues in order to clarify 
the universality of the FDEN model. This will allow us 
to clarify the applicability of our model to various tissue 
types. We believe that the FDEN model can represent 
other biological tissues consisting of a single tissue type, 
excluding non-soft tissues such as bone and tissue ex¬ 
hibiting plasticity, such as brain. Our previous nonlinear 
viscoelastic model with four parameters has al¬ 

ready been partially evaluated with breast tissues (fibro- 
glandur, fat, muscle) [13 HH]- We plan to evaluate the 
FDEN model using other tissues in future studies. 

In addition, this article did not utilize stress relax¬ 
ation and indentation tests, which are basic experiments 
to evaluate viscoelastic properties and stress-strain non¬ 
linearity, respectively. Parameter identification in the 
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TABLE I. Accuracy evaluation results for the present model. 


Equation and 
exp erimental data 

Equation ( 9a|)| 9b I and dynamic viscoelastic test 
E quation ( 13b I and creep results 
Equation! 16a I ( 16b|) and nonlinear measurement 


Sample number 
(Trial number) 

n 

64 (712) 

64 


Avg. of 
0.840 
0.997 
0.986 


Max. of 
0.902 
1.000 
1.000 


Min. oiR^ 
0.751 
0.950 
0.923 


S.D. of R^ 
0.008 
0.0073 
0.018 


TABLE 11. Fundamental statistics of the parameters. 


Sample 


Test 

Parameter 

number 

Avg. 

Max. 

Min. 

S.D. 

Dynamic viscoelastic test 

G 

11 

391.1 

518.2 

248.2 

111.2 

Dynamic viscoelastic test 

a 

11 

0.131 

0.146 

0.118 

0.011 

Nonlinearity measurement 

G 

64 

544.8 

1294 

341.8 

155.3 

Nonlinearity measurement 

G„ 

64 

8.547 

13.18 

5.26 

1.604 

Nonlinearity measurement 

Xb 

64 

0.121 

0.190 

0.076 

0.0022 

Nonlinearity measurement 

Gi 

64 

23.90 

39.64 

11.35 

6.206 


FDEN model using these tests are described in appen- Appendix A: Calculation of shear stress and strain 


dices 1^ and Stress relaxation analysis using fractional 
viscoelasticity was introduced in related studies [iQl |4T] , 
as well as in our work on human stretch applications |30j . 
Moreover, indentation (in the case of needle insertion 
and palapation for medical robotics) using the nonlinear 
model has been introduced in related studies [23H28] . 


Torque T applied to the sample, and the torsional an¬ 
gle 9 of the sample, were measured using a rheometer 
(AR-G2 or AR550; TA Instruments, New Castle, DE). 
From these measurements, the conventional shear strain 
X and conventional shear stress / were calculated using 
Eq. (Ala) and (Alb), respectively: 


V. CONCLUSION 

We proposed a simple empirical model using Frac¬ 
tional Dynamics and Exponential Nonlinearity (FDEN) 
to identify the rheological properties of soft biological 
tissue. The model is derived from detailed material mea¬ 
surements using samples isolated from porcine liver. We 
conducted dynamic viscoelastic tests and creep tests on 
liver samples using a rheometer. The experimental re¬ 
sults indicated that biological tissue has specific prop¬ 
erties, such as: i) power law increases in storage elastic 
modulus and loss elastic modulus with the same slope; ii) 
power law gain decrease and constant phase delay in the 
frequency domain over two decades; iii) log-log scale lin¬ 
earity between time and strain relationships under con¬ 
stant force; and iv) linear and log scale linearity between 
strain and stress relationships. Our simple FDEN model 
uses only three dependent parameters and represents the 
specific properties of soft biological tissue. 


(Alb) 

where d and R are the length and radius of the cylin¬ 
der (ref. Figl|), respectively. The mean stress and strain 
on the sample (half values of outer stress and strain on 
the sample) are referred to in the experimental results, 
because they are adequate for consideration of the non¬ 
linear properties. R was 20 [mm] and d was 5 [mm] in 
the experimental setup of this paper. 


Appendix B: Parameter dependency 


We modeled nonlinear properties of soft biological tis¬ 
sue bas ed on thes e resu lts and considerations, as shown 
in Eq. (Bla) and (Bib). 
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Gx = f {a: < Xh) 
= f {x > Xb} 


(Bla) 

(Bib) 


where x is strain and / is stress. G is linear stiffness, Xb 
is the boundary strain, G„ is nonlinear stiffness, Gi is the 
dependent parameter, and e is Napier’s constant. Each 
parameter should fulfill the condition that the exponen¬ 
tial curve (Bib is a tangent to the straight line (Bla) 
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at X = Xb- (B2a) and (B2b| is derived by derivation of 
equation (Bla) and (Bib); 


df_ 

dx 


^ = G {x < Xb} 
ax 

Gn{Gie^^^) {x>Xb} 


(B2a) 

(B2b) 


where G\ G” and w are variables; a and G are param¬ 
eters. 

We obtained the set of G’ and G” at each value for fre¬ 
quency u) from the experiment. We identified the param¬ 
eter from these data using EKF for system identification. 
The system identification in EKF is generally described 
as follows; 


The following equat ion m ay th us be fulhlled by plug¬ 
ging X = Xb into (|Bla| and (iBlb). 


= Gxb (B3) 

Moreover, the following equation may be fulfilled by 
plugging X = Xb into (B2a) and (B2b). 


= G 


(B4) 


By plugging the left side of equation (B3| into (B4), 


GjiXb — 1 


Then, 


Xb = 


G„ 


By plugging (B6| into (B4); 


(B5) 


(B6) 


9k+i = fiOkji’k) 
Vk = g{dk,Ck) 


(C2a) 

(C2b) 


where k = 0, 1, 2, represents the discrete iteration 
index (number of data set in this case), 0 is the n- 
dimensional state vector, ip is the n-dimensional system 
noise vector, y is the p-dimensional observation vector, 
C is the p-dimensional observation noise vector, and f() 
and g() are the nonlinear vector functions. In the theory 
of state-space, (C2a) and (C2b) are known as the sys¬ 


tem model (or state model) and the observation model, 
respectively. 

Parameter vector is regarded as a state vector in EKF 
for system identihcation. Where the state vector (param¬ 
eter vector) 6* is a constant vector and the observation 
noise vector ^ is a Gaussian white noise with zero mean, 


(C2a) and (C2b) are represented as: 


0k+l — dOk 
Vk = h[0k) + Ck 


(C3a) 

(C3b) 


G„G,e = G 


Then, 


G, = 


G 

GnC 


(B7) 


(B8) 


Each parameter may then fulhll the above relationship, 
particularly (B6) and ((B8|. The stress-strain relation¬ 
ship with several G„ in our model is described in Fig|lO| 
to represent the meaning of parameter constrain between 
G„ and Xb- 


Appendix C: Extended Kalman Filter (EKF) for 
Dynamic viscoelastic test 

This section shows the methodology used to identify 
the paramerter described in section [III A[ The model for 
dynamic viscoelastic test was as follows from equations 


where I is the identity matrix and h() is the nonlin¬ 
ear vector function. In the case of system identification 
for dynamic viscoelastic test, the state vector (parameter 
vector), 9 observation vector y and the nonlinear vector 


1000 



— S-S ctuve(Gn=15) 
Boundarypoint(Gn=15) 
S-S cuive(Gn=12.5) 
Boundary point(Gn=12.5) 

-S-S curve(Gn=10) 

Boundary point(Gn=10) 
S-S cuive(Gn=7.5) 
Boundary point(Gn=7.5) 
S-S curve(Gn=5) 
Boundary point(Gn=5) 


0.3 

Strain 


0.5 


0.6 


log G' = a logo; -|- log Gg. 
log G" = a log uj -\- log G'o 


(Cla) 

(Clb) 


FIG. 10. The stress-strain relationship with several Gn in our 
model to provide the meaning of parameter constrain between 
Gn and Xb- G was set to 1000 Pa. 
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function h() are particulary regarded as follows; 


a 

G 


y = 
h{9) = 


logG" 
logG" 

alogoj + log(Gcos(|Q;)) 
a logw + log(Gsin(|Q!)) 


(C4a) 

(C4b) 

(C4c) 


The EKF algorithm (ref. [BO]) using (C4a)-( |C4c ) was 
applied to identify the parameter from the data set. It 
was not necessary to set initial values for each parameter 
9q, meaning that 9o was zero vector. 


Appendix D: Extended Kalman Filter (EKF) for 
Nonlinearity measurement 

This section shows the methodology used to identify 
the paramerter described in section ElDi The model 
for nonlinearity measurment is as follows from equations 

(pal-p^; 


Gr, 


Gxc = fc {Xc < 7 ^}- (Dla) 

log(7^) = log(/c) {xc > 7^}. (Dlb) 

Lt7|,C Ltti 


where fc and Xc are variables, G and Gn are parame¬ 
ters, e is Neipias contant. We obtained the set of fc at 
each strain Xc from the experiment. We identified the 
parameter from these data using EKF for system identi- 
ficatrion. The algorithm to identify the parameter is the 
same as in Appendix particularly equations ( |C2a )- 
(C3b). In the case of system identihcation for nonlin¬ 


earity measurment, the state vector (parameter vector) 
0, observation vector y and nonlinear vector function h() 
are regarded as follows; 


used only low strain data (three data set from minimum 
strain) for approximation of the parameter G, while we 
used only high strain data (three data set from maxi¬ 
mum strain) for approximation of parameter G„. These 
approximations of the parameters can be identified using 
LSM -simple linear regression of equations (Dla) and 
dDlbl— . 


Appendix E: Relaxation test 

A model equation of stress in relaxation tests can be 
devised as follows. We assumed that equation @ is valid 
for a relaxation test, while nonlinearity was evaluated by 
a series of relaxation tests under several applied strains. 
Specihcally, equation ([^ becomes ( Ela[ ) if (|^ is solved 
for the conditions of the relaxation test. Here, the applied 
strain is constant Xc- Equation (Elbl is derived from the 
log-log transformation of (Ela). 


f = GT{l+r)xct-^{= fct-^. 
log / = -alogt-blog/c. 


(Ela) 

(Elb) 


where / is stress and t is time; Xc is constant strain; 
G is linear stiffness at each strain, and r() is the gamma 
function, fc is the coefficient determining the strain value 
as a parameter, which is defined as follows: 


fc = Gr(l -I- a)xc. 


(E2) 


The LSM algorithm can be used to identify the pa¬ 
rameters of equation (Elb) -linear regression- for each 


sample. We calc ulate d the other independent parameter 
G via equation ( |E2[ ). Nonlinear properties of samples 
can be investigated based on a series of relaxation tests 
under several applied strains. Specifically, nonlinearity 
measures the relationship between the constant applied 
strain Xc and the stress coefficient fc in the series of re¬ 
laxation tests under several strains. 


G 

G„ 


y = 


h{9) = 


f {x < 

ogf {x-^ 
Gx 




(D2a) 

cfcij 


(D2b) 


1-1 

V A 

(D2c) 


The EKF algorithm (ref. [BD]) using (D2a)-(D2c) was 


applied to identify the parameter from the data set. Each 
data set affected a single term of the vector, where the 
upper term for the vectors was updated via low strain 
data, and the lower term for the vectors was updated 
via high strain data. Initial values for each parameter 
00 needed to be explicitly set in the case of nonlinearity 
measurement. Therefore, we first approximated the pa¬ 
rameters to set initial values of each parameter 9q. We 


Appendix F: Indentation test 

A model equation for indentation test -the reac¬ 
tion force measurement during constant velocity strain 
change- can be theoretically calculated as follows. It 
should be noted that the sensitivity of parameter a from 
the steady-state of reaction force / is generally low in ex¬ 
periments with soft biological tissue. Linear stiffness G 
and nonlinear stiffness G„ can thus be identified from re¬ 
action force on indentation test, when the value of param¬ 
eter a is roughly known. We assumed that equations (la) 
and © -the FDEN model- are valid for indentation 
test, when the value of parameter a is already known. We 
collected time seies data for force f(t) on the conditions 
of indentation test. Here, the applied strain is x(t)= Vot. 
Equations (@- @ become identical to static problems 












































17 


such as (16a)-(16bl and (Fla)-(Flal when the fractional 


integrated stress /’ is considered as follows: 


Gx = f = D^-^^f {x<x4.(Fla) 
GnX + logGi = log/' = D^~°'Hogf{xc > a;h}.(Flb) 


where D(a) refers to ath-order derivative and D(-q:) 
refers to ath-order integral. The numerical fractional in¬ 
tegration, which is necessary in the above calculation is 
introduced in various studies (for example, m). Param¬ 
eters (G, G„, Gi and Xb) of equations (Flal and (Fib) 
can be identified via the same method introduced in Sec¬ 
tion 3.C, while we use fractional integrated stress /’ on 
behalf of stress /. 
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